arXiv: 1505.06616vl [hep-lat] 25 May 2015 


Prepared for submission to JHEP 


The Charmonium Potential at Non-Zero 
Temperature 


Chris Allton,“ Wynne Evans , a,b Pietro Giudice, c and Jon-lvar Skullerud d 

“ Swansea University, Swansea, United Kingdom 
b Bern University, Bern, Switzerland 
c Universitat Munster, Munster, Germany 
d University of Maynooth, Maynooth, Ireland 

E-mail: pyevans@swansea.ac.uk, c.allton@swansea.ac.uk, 
p.giudice@uni-muenster.de, jonivar@thphys.nuim.ie 

ABSTRACT: The potential between charm and anti-charm quarks is calculated non-perturb- 
atively using physical, rather than static quarks at temperatures on both sides of the de¬ 
confinement transition Tq, using a lattice simulation with 2+1 dynamical quark flavours. 
We used the HAL QCD time-dependent method, originally developed for inter-nucleon po¬ 
tentials. Our lattices are anisotropic, with temporal lattice spacing less than the spatial 
one which enhances the information content of our correlators at each temperature. Local- 
extended charmonium correlators were calculated efficiently by contracting propagators in 
momentum rather than coordinate space. We find no significant variation in the central 
potential for temperatures in the confined phase. As the temperature increases into the 
deconfinement phase, the potential flattens, consistent with the expected weakening inter¬ 
action. We fit the potential to both the (a) Cornell and (b) Debye-screened potential forms, 
with the latter better reproducing the data. The zero temperature string tension obtained 
from (a) agrees with results obtained elsewhere, and it decreases with temperature, but at a 
slower rate than from the static quark approximation. The Debye mass from (b) is close to 
zero for small temperatures, but starts to increase rapidly around T c - The spin-dependent 
potential is found to have a repulsive core and a distinct temperature dependence above Tq 
at distances 1 fm. 
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1 Introduction 

Shortly after the discovery of the particle, it was realised that charmonium bound 
states could be very well described using the Schrodinger equation with a phenomenological 
potential between the quark and antiquark [1-3]. This was one of the factors leading to QCD 
being accepted as the theory of the strong interaction, with charmonium as the ‘hydrogen 
atom’ of QCD. The validity of a potential model has since been formally established using 
an effective theory — potential nonrelativistic QCD (pNRQCD) — which is obtained by 
integrating out degrees of freedom with momentum above the typical binding energy of the 
system [4, 5]. The resulting potential for infinitely heavy (static) quarks can be shown to 
be equivalent to the one extracted from the Wilson loop, which has been computed non- 
perturbatively from first principles on the lattice. Recently, the potential between quarks 
with finite mass has also been computed from lattice QCD by ‘reverse engineering’ the 
Schrodinger equation [6-11] using the HAL QCD method [12]. 

At non-zero temperature, a potential model incorporating colour-Debye screening has 
been used to predict the dissociation of charmonium states as a signal for the formation 
of quark-gluon plasma (QGP) [13]. Since that seminal paper, a substantial experimental 
effort has been invested in the study of ,//'(/> suppression at SPS, RHIC and the LHC [14— 
17], and a number of studies have been carried out using potential models for charmonium 
systems at high temperature [18-20]. However, unlike at zero temperature, there was no 
rigorous proof of the validity of potential models for static quarks and hence no agreement 
on what to use for the inter quark potential. Different groups have used, for example, the 
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free energy [21, 22] or the internal energy [23] of static quarks as computed on the lattice, 
or a combination of these [24, 25]. 

Recently, a series of effective theories has been developed, depending on hierarchies 
such as 


M q » T > g 2 M q > gT > g 4 M q or M q » g 2 M q > T > gT > g*M q , (1.1) 

where M q is the heavy quark mass, T is the temperature and g is the gauge coupling. A 
common feature of these theories is the appearance of an imaginary part of the potential, 
resulting from Landau damping [26, 27]. In certain parameter ranges this term can be more 
important for charmonium suppression than the Debye screening encoded in the real part. 

As in the zero temperature case, non-perturbative (lattice) calculations of the Wilson 
loop have been used to extract the static inter-quark potential at non-zero temperature 
[28-31]. 

In this work we further extend these non-perturbative lattice calculations to the finite 
quark mass case. We will assume the validity of the Schrodinger equation and potential 
description for charm quarks, and derive the charmonium potential at non-zero temper¬ 
ature directly from charmonium correlators following the HAL QCD method. We restrict 
ourselves to the real part of the potential and will be following the ‘time-dependent’ method 
introduced in [32], This method was used in [33] but is distinct from the ‘fitting’ method 
used in [34]. In the fitting method, local-extended correlators are first fitted to exponen¬ 
tials at large Euclidean time, r, to extract the Nambu-Bethe-Salpeter (NBS) ground state 
wave function. The NBS wave function is then used, in conjunction with the Schrodinger 
equation, to reverse-engineer the potential. The fitting method is well understood from 
a theoretical point of view since it relies on conventional fitting techniques. However, at 
non-zero temperature, where the temporal range of the correlators is limited, it suffers from 
familiar limitations: higher excited states still contribute to the correlator at the largest 
available r, making fits unreliable. 

In [35] the first calculation of the charmonium potential at finite mass and non-zero 
temperature was performed, but this used the quenched approximation. We performed 
thermal studies of the charmonium potential using our two-flavour, 1st generation FASTSUM 
ensembles in our previous study [33, 34], In the work presented here, we extend this by using 
our 2nd generation ensembles which have 2+1 flavours with finer, larger lattices and light 
quarks closer to their physical masses. We find that the potential does not vary significantly 
for temperatures below the crossover temperature, Tc, but that it clearly flattens above 
T(j. Using the Cornell form of the potential, we determine the temperature dependent 
string tension which decreases as T increases as expected. We also fit our data to a Debye- 
screened potential form to determine the Debye mass. The spin-dependent potential is also 
determined and we find thermal effects for T 'ft Tq . An early version of this work appears 
in [36], 

In section 2, the HAL QCD time-dependent method is reviewed, and we discuss the 
effect of the backward moving states. Our simulation details are outlined in section 3 and 
section 4 presents our results. 
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Figure 1. Quark propagator diagram of a local-extended correlator. The quark and anti-quark 
are created at the source with no separation and annihilated at the sink with a separation r. 


2 The Method 


The HAL QCD time-dependent method of extracting the potential from the lattice [32] takes 
local-extended correlators as input, see Figure 1. These are constructed from creation and 
annihilation operators that have the form 

Jr(x; r) = q(x ) T U(x, x+r) g(x+r), (2-1) 


where r is the separation between the quark and anti-quark fields q and q, x is the space- 
time point (x, r) and T is a monomial of gamma matrices used to generate pseudoscalar 
( rj c ) channels or vector («//?/>). [7(x,x + r) is the gauge connection between x and x + r 
required for gauge invariance. 

The local-extended charmonium correlator is 


Cr(r, r) = J r (x, r; r) 4(0; 0)), (2.2) 

X 


where the sum over the spatial coordinate at the sink, x, projects the momentum of the 
state to zero. 

The local-extended correlator can also be expressed as a sum over the eigenstates of 
the Hamiltonian with eigenvalues Ej, 


Cr(r,r) 


^ 2 Ei 

3 J 



e ~Ej(N T 



(2.3) 


where the sum is over the states j with the same Lorentz transformation properties as the 
operator Jp, N T is the number of lattice points in the temporal direction and Vh'( r ) is the 
charmonium wave function. We now consider only the forward-moving contribution to the 
correlator (the effect of ignoring the backward-moving contribution is discussed later): 


Cr(r,r) 


E 


^(0)^(r) Ei r 
2 Ej 


E V 4(r)e EjT , 


(2.4) 


where we have defined the unnormalised wavefunction, \E r J (r) = 'il>*(0)'ipj(r)/2Ej. We treat 
the charm quark non-relativistically due to its large mass, and assume \F,(r) obeys the 
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Schrodinger equation, 


(2.5) 


(-^ + Mr)) * j (r) = E j * j (r), 

where Vy(r) is the desired potential for the channel T, n = m c /2 is the reduced mass, and 
we only consider S-wave states. Taking the time derivative of (2.4) and using (2.5), we 
obtain, 


<9Cr(r,r) 

8 t 


3 


- yr ( r )) ^( r ) e 

(5“^ r(r) ) Cr(r,T) 


( 2 . 6 ) 


which can be trivially rearranged to yield the potential, 

Vi ' (r) = crib)(l^) c ’ r<r ' T) 


(2.7) 


We highlight the fact that Vy(r) from (2.7) has an implicit r dependence which must be 
averaged over, see section 4. 

On the lattice, the Laplacian in (2.7) is approximated as follows, 


V?/(r) 


' d 2 2d' 

dir 2 r dr 


fir) 



,'r+a s 


25 r ' r + 3r',r—a 3 $r' ,r+a s &r',r—a s 


+ 


ra. 


fir'), 

( 2 . 8 ) 


where we have relied on the approximate rotational symmetry of the lattice. 

The time derivative in (2.7) can be approximated by the naive finite temporal difference, 


d n \ , \fi T + a r) ~ fir - a T ) 

-* - 2^- 


(2-9) 


but, as we will see in section 4.1, this approximation is particularly poor near the temporal 
centre of the lattice because of contamination by the backward mover which we neglected 
in the above derivation going from (2.3) to (2.4). For large temperatures, corresponding to 
lattices with a small temporal extent, this is especially problematic because the uncontam¬ 
inated region may become vanishingly small. 

Fortunately the expression [37] 


_ 1 w ( CyjT - a T ) + y/Cr(r - a T ) 2 - C T (N T a r /2y \ 

° T 2 ° g \ Ct (t + or) + yf Cy{r + a T ) 2 - C v iN T a T /2y) ’ 


( 2 . 10 ) 


recovers the exact ground state energy, Eq, even in the case where there is a backward 
mover as in (2.3) (in the absence of excited states). 

In the next section, we test the improved temporal difference based on (2.10), 


d_ 

df 


C'r(r) 
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N s 

N t 

T(MeV) 

T/T c 

A^cfg 

24 

128 

44 

0.24 

250 

24 

40 

141 

0.76 

500 

24 

36 

156 

0.84 

500 

24 

32 

176 

0.95 

1000 

24 

28 

201 

1.09 

1000 

24 

24 

235 

1.27 

1000 

24 

20 

281 

1.52 

1000 

24 

16 

352 

1.90 

1000 

32 

32 

176 

0.95 

500 

32 

24 

235 

1.27 

500 


Table 1. The spatial and temporal extents, N s and N T , the temperature in units of MeV and Tc, 
and the number of configurations, N c { g , of the ensembles used. Note that for two temperatures, we 
have two different spatial volumes to study finite volume effects. 


to define the potential, (2.7). We find it has significantly reduced contamination from the 
backward mover compared to the naive expression (2.9). 

Once the potential for both the pseudoscalar (PS) and vector (V) channels have been 
determined from the method described above, the central and spin-dependent potentials, 
Vq and Vg, can be derived as follows. The leading order terms in the velocity expansion of 
the interquark potential for S-wave states can be expressed as [38, 39], 

^r(r) = Vc(r) + y s (r)ai-s 2 , (2.12) 

where * 1.2 are the quark spins. Using si • S 2 = —3/4, 1/4 for the PS and V channels 
respectively, Vc(r) and Vg(r) can be obtained from the Vps,V potentials using the following 
expressions, 

Veto = Jups + ^Uy, U s (r) = U v - ^PS- (2.13) 

3 Simulation Details 

We use our FASTSUM Collaboration’s 2nd generation configurations for this analysis [40]. A 
Symanzik gauge action and a stout-smeared clover fermion action were used with parameters 
set by the Hadron Spectrum Collaboration (hsc) in [41, 42], There are 2+1 flavours of dy¬ 
namical quarks with the two (degenerate) light flavours corresponding to AI n = 392(4) MeV 
and the third dynamical quark set to the strange quark mass. Anisotropic lattices were 
employed with an anisotropy of £ = a s /a T ~ 3.5, with a s — 0.123 fm and a T 1 ~ 5.63 GeV. 
The lattice spacings are fixed and we vary the temperature, T = (a r V T ) -1 , by adjusting the 
number of temporal points, N T . The ensemble parameters are listed in Table 3 showing that 
we study temperatures in both the confined and deconfined phases. The “zero” temperature 
(i.e. N r = 128) ensemble was kindly provided by the HSC Collaboration. Our main spatial 
volume is (3fm) 3 , corresponding to N s = 24, and we have two temperatures with a (4fm) 3 
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volume ( N s = 32) to enable us to check finite volume effects. The deconfinement crossover 
temperature, Ty, was determined from the inflection point of the Polyakov loop [40, 43]. 

The (non-dynamical) charm quark was calculated with the same (relativistic) action 
used for the three light dynamical quarks with its mass set by tuning the PS state to the 
experimental rj c mass while simultaneously maintaining the anisotropy [44], As in [7], the 
quark mass is defined as My/ 2 where My is the mass of the charmonium vector channel 
ground state. Hence the reduced mass, /i = My/ 4, see (2.5). 

We chose to gauge fix our configurations to the Coulomb gauge, and then replace the 
gauge connection, U(x,x + r), in (2.1) by unity. We used the highly optimised Fourier- 
accelerated gauge fixing procedure of [45]. The coordinate space quark propagators were 
calculated using the Chroma software suite [46] and then tied together in momentum space 
using a bespoke program to obtain correlators more efficiently, see the Appendix for details. 

4 Results 

4.1 Central Potential 

Local-extended correlators, (2.2), corresponding to on-axis quark separations were gener¬ 
ated from the ensembles outlined in section 3. Correlators corresponding to quark separa¬ 
tions of the same magnitude, |r|, were averaged giving 13 unique separations for N s = 24. 
The set of PS correlators for the 0.76Tc ensemble is shown in Figure 2 for all available r. 
As one would expect, the signal decreases as the quark separation, r, increases. Assuming 
ground state dominance of the correlator, this follows from the monotonic property of the 
S-wave wavefunction. 

The N s = 32 ensembles, listed in Table 3, provide a means to investigate the volume 
dependence of the correlators, and hence also that of the potentials. Figure 3 shows the 
ratio of N s = 24 to N s = 32 local-extended correlators for the 1.27T C (i.e. the N T = 24) 
case. From Figure 3, the correlator has no volume dependence for 0 < r/a s < 10, the 
r/a s = 11 case shows some effects and the r/a s = 12 is clearly highly sensitive to finite 
volume effects. Consequently, (due to the nearest-neighbour representation of the Laplacian 
in (2.8)) we report the potential for r/a s < 9, i.e. r < l.lfm, in the following where it is 
free from finite volume effects. 

We note that the lattice version of the Laplacian in (2.8) has greatest discretisation 
error near the origin. Hence we do not include the r/a s = 1 point in our fits to the potential 
in section 4.2. 

In Figure 4 the PS potential is plotted for the 0.76Tc correlators using the naive 
time derivative (2.9) in (2.7). We note that in the HAL QCD method, the potential has 
a r pseudo-dependence which will be averaged over. Figures 5 and 6 show the spatial 
and temporal contributions of (2.7) respectively. As can been seen, at large r the spatial 
derivative contribution is stable but the temporal derivative contribution increases near 
the mid-point of the lattice, r ~ N T / 2, which produces a corresponding decrease in the 
potential at these r values in Figure 4. In Figure 7 we use the improved temporal term 
(2.10) which clearly resolves the issue for the pseudoscalar case, implying that it was caused 
by the backward mover as discussed in section 2. We have checked that the same is true in 
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Figure 2. Local-extended charmonium correlators for the PS channel for all possible on-axis 
separations of the 0.76T C (N T = 40) ensemble. Error bars are smaller than the symbols. 
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Figure 3. The ratio of the N s = 24 to N s = 32 local-extended correlation functions for 1.277), (i.e. 
N t = 24). In the lower pane, all separations 0 < |r/o T | < 12 are shown. The upper pane shows a 
closeup for 0 < |r/o T | < 10 with the points shifted horizontally for clarity. 
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Figure 4. The potential for the 0.76Xc pseudoscalar channel plotted as a function of r for each 
quark separation |r|, using (2.7) with the naive time derivative (2.9). 


the vector channel. We also note that as r —> oo there is the expected convergence towards 
the (negative) value of the PS mass r/ c = 2.9804(l)GeV [47]. We note that in the region 
where (2.7) is valid, i.e. in the absence of a backward mover, Eq is equivalent to dCr/dr 
only in the limit of negligible excited state contribution. However, noting the very similar 
small-r behaviour in Figures 6 and 7 (where the excited states’ contributions are largest), 
it is clear that the Eq definition works in practice even in the presence of excited states. 
For these reasons, we will always use the improved temporal form (2.10) in the following. 

The central potential is determined by combining the pseudoscalar and vector potentials 
according to (2.13) and is shown in Figure 8 for all temperatures. From these figures, it is 
clear that the 1.52Tc and 1.90Tc data do not stabilise. Hence, we include only the data 
up to temperatures of 1.27Tc in the following. 

We remove the r pseudo-dependence from the potentials by performing a correlated 
fit to a constant using the r ranges shown in Table 2. Two time ranges are chosen to 
elucidate the systematic error from the choice of fit range. The “best” range is chosen to be 
the best available for each temperature, and the “lower” range is chosen to match the best 
range of the next higher temperature. In this way, direct comparisons can be made between 
neighbouring temperatures. Using these time ranges, we obtain the temperature-dependent 
potential shown in Figure 9. The circles correspond to the best range and triangles the lower 
range. To aid this comparison between temperatures, we include five upper insert plots in 
Figure 9 which show a closeup of the data at every second separation, i.e. r = a s , 3 a s ,... 9a s . 
The vertical range of all the insert plots is identical (0.6 GeV). We emphasize from the 
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Figure 5. The spatial contribution, (l/C'r)(V^C'r/2/r), to the potential (see (2.7)) for the 0.76Tc 
pseudoscalar channel, plotted as a function of r for each quark separation |r|. 
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Figure 6. The temporal contribution, (l/Cr)(9Cr/9r), to the potential (see (2.7)), using the 
naive form, (2.9), for the 0.76Tc pseudoscalar channel, plotted as a function of r for each quark 
separation |r|. 
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Figure 7. The temporal contribution, (l/Cr)(c?Cr/c?T), to the potential, using the improved 
temporal form, Eq, see (2.10), for the 0.76Tc pseudoscalar channel for each quark separation |r|. 


T/T c 

N r 

Best 

range 

Lower 

range 

0.24 

128 

30 

- 63 

15 - 

- 19 

0.76 

40 

15 

- 19 

12 - 

- 17 

0.84 

36 

12 

- 17 

11 - 

- 15 

0.95 

32 

11 

- 15 

11 - 

- 13 

1.09 

28 

11 

- 13 

9 - 

11 

1.27 

24 

9 - 

- 11 

N/ 

'A 


Table 2. Fitting ranges used to obtain the r—independent potentials. The “best” range covers the 
plateau region for each temperature and gives our best fit. The “lower” range is the same as the 
best range for the next higher temperature, allowing a direct comparison between temperatures. 


discussion above, a direct comparison between neighbouring temperatures can be made by 
comparing the triangles at one temperature with the circles at the next higher temperature. 
We therefore conclude that the T ^ Tq potentials have no significant temperature effects 
that our data can discern for any separation, but that there is a significant flattening of the 
potential at moderate to large separations, r ^ 0.3 fm, for T ^Tq. 

In Table 3, the central potential values are listed with the systematic (from the t fitting 
procedure described above) and statistical errors combined additively into a single error bar. 
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Figure 8. The central potential for all temperatures studied as a function of r for each quark 
separation, |r|. 


11 














r [fm] 
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Figure 9. The central potential for various temperatures as a function of separation. Two different 
r ranges were used to estimate the systematic error for all but the largest temperature, see Table 
2. The circles represent the “best” fit range in r for each temperature, and the triangles the “lower” 
fit range (which is the next highest temperature’s best range). This allows a direct comparison at 
different temperatures to be made. The five upper graphs are closeups of the separations r/a s = 
1, 3,... 9 with a common vertical range of 0.6 GeV. All points are shifted slightly horizontally for 
clarity. 


4.2 The Cornell Potential, String Tension and Debye Screening 

In Figure 10 the central potential is plotted with the combined statistical and systematic 
errors as listed in Table 3. For comparison, the Cornell potential [1-3], 

V(r, T) = + a (T)r + C, (4.1) 

r 

is also shown using continuum parameters a c = 7r/12 [48] and \J~o = 0.445 GeV as used 
in [49] (with C adjusted to overlie our data). It was shown in [49] that these parameters 
reproduce the properties of the lowest lying states in both charmonium and bottomonium 
very accurately. Since our “zero” temperature (i.e. T = 0.24Tc) data follow this established 
Cornell potential extremely well, they provide strong evidence that our method is extracting 
the correct physics. 

The temperature dependence of the central potential is further studied by fitting the 
Vq data to the Cornell function (4.1) for both the “best and “lower” r ranges, see Table 
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r [fm] 

0.24 T c 

0.76 T c 

[GeV] 
0.84 T c 

0.95 T c 

1.09 T c 

1.27 T c 

0.123 

2.58(6) 

2.63(3) 

2.66(2) 

2.67(2) 

2.67(5) 

2.703(2) 

0.246 

2.87(6) 

2.92(4) 

2.95(2) 

2.96(2) 

2.94(6) 

2.967(5) 

0.369 

3.06(6) 

3.12(5) 

3.14(1) 

3.15(3) 

3.11(7) 

3.121(7) 

0.492 

3.23(5) 

3.28(5) 

3.29(1) 

3.29(3) 

3.24(6) 

3.23(1) 

0.615 

3.38(3) 

3.43(5) 

3.42(1) 

3.41(4) 

3.35(4) 

3.30(1) 

0.738 

3.51(1) 

3.55(5) 

3.54(2) 

3.50(4) 

3.43(2) 

3.35(2) 

0.861 

3.63(1) 

3.66(5) 

3.63(3) 

3.58(4) 

3.50(2) 

3.38(2) 

0.984 

3.73(3) 

3.75(5) 

3.72(5) 

3.64(4) 

3.57(4) 

3.40(3) 

1.107 

3.82(3) 

3.85(5) 

3.80(6) 

3.70(4) 

3.64(8) 

3.41(3) 


Table 3. The central potential data as shown in Figure 10. The statistical uncertainty has been 
combined additively with the systematic error from the r—range choice as described in the text, 
but note that the errors for the 1.27Tc data are statistical only. 


2. We restrict the fit to separations in the range 2 a s < r < 9 a s due to the systematics at 
very small and large separations discussed in section 4.1. The parameters from these fits 
are shown in Table 4 and the resultant string tensions are plotted in Figure 11. We also 
show the fits in Figure 10 for the coldest and hottest temperatures. 

At face value, there is a clear temperature dependence in the string tension as displayed 
in Figure 11. However, using the strict criteria outlined in section 4.1, the neighbouring tem¬ 
perature’s “best” and “lower” string tensions overlap within errors, and so higher statistics 
are required to properly decouple this thermal effect from the systematics in this quantity. 

A calculation of the string tension at zero temperature was performed in [9] with 
nearly physical light dynamical quarks using the HAL QCD method, where they found a = 
394(7) MeV. In studies using the static quark potential, the string tension displays clear 
temperature effects well below Tc, see e.g. [50] (in contrast to Figure 11). However [50] 
study pure gauge SU (3) where the confining string tension is essentially an order parameter 
which is therefore zero for T > Tq,. 
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Figure 10. The central potential with combined statistical and systematic errors for various 
temperatures (noting that the 1.277c data has statistical errors only). Data points are shifted 
slightly horizontally for clarity. Examples of the fits to our data are shown using the Cornell (4.1) 
and Debye-screened (4.2) potentials. These used the restricted range 2 a s < r < 9 a s as described 
in the text. The Cornell form from (4.1) with the continuum parameters a c = 7t/12 [48] and 
yfo = 0.445 GeV [49] is shown which agrees extremely well with the Cornell fit to our 0.24T C data. 


In [51], an alternative to the Cornell potential (4.1) was proposed for finite temperature 
where there is colour-Debye screening of the colour sources. This has the form, 

V(r,T) = - as ^ T l e -^D(T)r + ^ = 0) / 1 _ e -m D (T)r\ + ^ ( 4 . 2 ) 

r mD\T) V / 

where m/j(T) is the Debye screening mass. This functional form has the feature that V(r, T ) 
remains finite as r —> oo. 

We have performed fits of the central potential for each temperature using (4.2) with 
three fitting parameters, a, mo and C. We fixed a to its “zero” temperature value (obtained 
with our T = 0.24Tc ensemble) of 434 MeV, see Table 4. As the temperature increases, the 
screened form, (4.2), fits the data better than the Cornell form, (4.1). This can be seen in 
Figure 10 where both the screened and Cornell fit for the hottest temperature are shown. 
For the coldest temperature, the screened and Cornell fits are indistinguishable. 

The results of these fits are shown in Table 4 and the resultant Debye masses are 
plotted in Figure 12. At low temperatures, mo ~ 0 and it then has a rapid increase 
around Tq, in agreement with expectations. In [52], the Debye mass was calculated from 
a two-flavour lattice calculation of the static quark-antiquark free energy. They found very 
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T/T c 

r-range 

Cornell Fit (4.1) 

Screened Fit (4.2) 



\fo 

[MeV] 

«c 

m D 

[MeV] 

«s 

0.24 

0.24 

30-63 

15-19 

434(7) 

390(20) 

0.31(2) 

0.39(5) 

3(8) 

50(30) 

0.311(9) 

0.33(3) 

0.76 

0.76 

15-19 

12-17 

410(20) 

405(14) 

0.36(3) 

0.41(3) 

30(20) 

40(20) 

0.34(2) 

0.38(2) 

0.84 

0.84 

12-17 

11-15 

378(13) 

340(20) 

0.41(3) 

0.46(3) 

70(20) 

110(20) 

0.35(2) 

0.36(2) 

0.95 

0.95 

11-15 

11-13 

331(11) 

340(20) 

0.48(2) 

0.49(3) 

130(20) 

120(30) 

0.37(2) 

0.39(2) 

1.09 

1.09 

11-13 

9-11 

320(20) 

220(20) 

0.43(3) 

0.58(3) 

140(20) 

270(30) 

0.31(2) 

0.42(3) 

1.27 

9-11 

180(30) 

0.53(3) 

340(40) 

0.37(3) 


Table 4. The results of fitting the central potentials to both the Cornell form, (4.1), and the Debye- 
screened form, (4.2), as described in the text. Note that, following the procedure outlined in section 
4.1, two t— ranges were chosen to allow a direct comparison between neighbouring temperatures. 


similar behaviour to Figure 12, with mp increasing rapidly around Tq to around 400 MeV 
at 1.2 Tq. 

4.3 Spin Dependent Potential 

We also obtained the (r—dependent) spin-dependent potentials for the different tempera¬ 
tures by combining the pseudoscalar and vector time-slice potentials from each ensemble 
according to (2.13). The r—dependence was removed using the same procedure as in the 
central potential case to obtain Figure 13. 

Taken at face value, we see a strongly repulsive core, but this has to be qualified by the 
systematics in the r/a s = 1 data point as discussed in section 4.1. However, since the spin- 
dependent potential is the difference, Vy — Vps, the systematics at r = a s may cancel to 
some extent. Also note that [8-11] have found a repulsive core for this potential, but with a 
quenched calculation. Modelling the interaction via one-gluon exchange, the spin-dependent 
potential is a <5—function at the origin, so given the body of lattice results including this 
work, a finite-width, repulsive core appears to be the correct, non-perturbative result. Note 
also the work of [53] where a finite-width repulsive potential was obtained by including the 
running of the coupling in one-gluon exchange. 

For confined temperatures, T ^ Tp, the potential is clearly fiat for moderate to large 
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Figure 11. The string tension, cr, as a function of temperature obtained by fitting the potential 
data to the Cornell potential (4.1). Two different r—ranges “best” and “lower” were also used so 
that temperature effects can be uncovered, see Table 2. The points are shifted slightly horizontally 
for clarity. 


distances with no significant temperature dependence, see Figure 13. In our calculation, as 
in [8, 10] the asymptotic value for this confined phase potential is negative. However, when 
the spin-dependent potential is used to define dynamically the reduced mass, //, (see (2.5)) 
this potential tends to zero at large distances [9, 11] by definition. 

There is a clear temperature effect once the deconhned phase is reached with a distinct 
minimum at intermediate distances r ~ 0.4 fnr and significantly larger potential values at 
large distances r ^ 0.7 fnr compared to the same potential in the confined phase. 

4.4 Comparison with other methods 

In Figure 14 our central potentials from this work (i.e. using our Nf = 2 + 1 2nd generation 
ensembles) are compared with those obtained from our earlier, Nf = 2 1st generation sim¬ 
ulations [33]. The Nf = 2 potentials were also obtained with the HAL QCD time-dependent 
method, and have been shifted vertically in Figure 14 so that their r/a s = 1 data points 
coincide with that of the 0.24Xc potential. It is encouraging that the potential data points 
interpolate each other at small r, especially given that the lattice parameters and actions 
used in each simulation are quite different. For a given temperature the Nf = 2 + 1 central 
potentials are flatter at large r than those from the Nf = 2 simulation. This could be due 
to the inclusion of an extra sea quark that has the ability to screen the strong force between 
quarks, but further studies would be required to confirm this. 
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Figure 12. The Debye mass, mo, as a function of temperature obtained by fitting the potential 
data to the screened potential (4.2). Two different r—ranges “best” and “lower” were also used so 
that temperature effects can be uncovered, see Table 2. The points are shifted slightly horizontally 
for clarity. 


In Figure 15 we compare the central potentials with the static quark potentials calcu¬ 
lated from the Wilson lines [54] which were also obtained with 2+1 flavours. The static 
quark potential curves in Figure 15 are shown at the temperatures closest to those in this 
work, and have been shifted vertically so that their form can be compared to our result. 1 
While the higher temperature results agree fairly well between the two approaches, the 
lower temperature static data are steeper than our results. Further study would be required 
to determine if this difference is due to [54] using the infinite quark mass approximation 
rather than the physical charm quark, or to other systematic differences between the two 
approaches. 

5 Conclusions 

There is a significant body of theoretical work studying the interquark potential at non-zero 
temperature using model, perturbative and lattice (non-perturbative) approaches. However, 
until [33, 35], these lattice studies all used the static quark limit. This work improves upon 
these static calculations by considering quarks with finite mass, and thus represents a first- 
principles calculation of the charmonium potential of QCD at non-zero temperature. The 

1 This is justified since in the static limit the quark mass, which sets the overall scale, has been removed 
from the free energy calculation. 
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Figure 13. The spin-dependent potential with combined statistical and systematic errors for 
various temperatures. The points are shifted slightly horizontally for clarity. 

method we used is based on the HAL QCD time-dependent approach which obtains the 
potential directly from local-extended correlators. 

We do not observe any significant temperature dependence of the central potential 
below Tq, while there is a significant flattening above Tq, consistent with the expectation 
that the potential becomes deconfining. The string tension is calculated and we find a 
slower variation of this quantity with temperature than that found using the static quark 
approximation. Using the Debye-screened form for the potential (which fits our data better 
than the Cornell form), we determine the Debye mass which is found to be very small at 
low temperatures and then increase rapidly around Tq. This is, as far as we know, the 
first non-perturbative calculation of the Debye mass in charmonium. In the case of the 
spin-dependent potential, we similarly find no thermal modification for T ^ Tq, but a clear 
variation with temperature at large distances in the deconfined phase and evidence for a 
repulsive core. 

This work improves upon our earlier work [33, 34] by including a dynamical strange 
quark, and using lattices which are finer, with a larger volume, and have lighter, more 
physical u, d quarks. 
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Figure 14. A comparison of the central potentials obtained from this Nf = 2 + 1 work and our 
earlier Nf = 2 simulations [33]. The Nf = 2 data have two error bars, statistical (left) and system¬ 
atic (right), and are shifted vertically so that the potential at the first separation agrees with the 
Nf = 2 + 1, 0.24Tc values. All points have been shifted slightly horizontally for clarity. 
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A Appendix: Momentum Space Propagators 


Local-extended correlators can be obtained more efficiently by working with the quark 
propagators in momentum rather than coordinate space [55]. While this method (which is 
not our own) has been known for some time and used in many papers which have calculated 
potentials from the HAL QCD method and in studies of multi-baryon states [56], it does not 
appear to be in any publication, so we outline it here for reference. 

For a meson, the local-extended correlator, (2.2), in the gauge fixed case (with U(x, a:+r) 
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Figure 15. A comparison of the central potential obtained from this work with the static quark 
potential for comparable temperatures obtained in [54]. The static potential data have been shifted 
vertically for the sake of comparison. 


set to unity) can be written, 

CY (r, r) = — ^^Tr H - 1 (x+r,r: 0 , 0 )rZl _ 1 ( 0 , 0 :x,r)rT 

X 

= — ^Tr H~ 1 (x+r, r:0,0)r7 5 Zl~ 1 (x, r:0,0)75r^ 


(A.l) 


where D 1 (sink: source) is the quark propagator. This correlator can be written in terms 
of the Fourier transform of the quark propagators, 

D-\ y, r: 0 , 0 ) = 1 £ ^(q)^, (A. 2 ) 

q 


giving 


Cr(r,r) 


1 

V 


E Tr 


n '(pir-.r.D '(-pi-.-.r 1 


±£C(p,T)e‘r 

P 


e ipr 


(A. 3 ) 


where C(p,r) is the Fourier transform of the correlator, Cr(r,r). 

This implies that once Cr(p, r) is obtained, then the desired correlator, Cr(r, r), can be 
simply computed for any value of r using the single sum in (A. 3 ). This is computationally 
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more efficient than having to perform the trace in (A.l) for each x value before then also 
performing the sum over x. 
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